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Glossary 

This glossary consists of a list of terms used inside the paper: In mathematics, 
in probability, in engineering, and on occasion in physics. To clarify the seemingly 
confusing use of up to four different names for the same idea or concept, we have 
further added informal explanations spelling out the reasons behind the differences 
in current terminology from neighboring fields. 

Disclaimer: This glossary has the structure of four columns. A number of terms 
are listed line by line, and each line is followed by explanation. Some "terms" have 
up to four separate (yet commonly accepted) names. 

MATHEMATICS PROBABILITY ENGINEERING PHYSICS 

function random variable signal state 

(measurable) 

Mathematically, functions may map between any two sets, say, 
from X to Y; but if A is a probability space (typically called fi). 



Work supported in part by the U.S. National Science Foundation. Tlie full version with figures 
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MATHEMATICS PROBABILITY ENGINEERING PHYSICS 

it comes with a cr-algcbra B of measurable sets, and probability 
measure P. Elements E in B are called events, and P(E) 

the probability of E. Corresponding measurable functions with 
values in a vector space arc called random variables, a terminology 
which suggests a stochastic viewpoint. The function values of a 
random variable may represent the outcomes of an experiment, 
for example "throwing of a die." 

Yet, function theory is widely used in engineering where 
functions are typically thought of as signal. In this case, X may 
be the real line for time, or M.'^. Engineers visualize functions as 
signals. A particular signal may have a stochastic component, 
and this feature simply introduces an extra stochastic variable 
into the "signal," for example noise. 

Turning to physics, in our present application, the physical 
functions will be typically be in some L^-space, and i^-functions 
with unit norm represent quantum mechanical "states." 

sequence (inch random walk time-series measurement 
vector- valued) 

Mathematically, a sequence is a function definc;(l on the integers 
Z or on subsets of Z, for example the natural numbers N. Hence, 
if time is discrete, this to the engineer represents a time series, 
such as a speech signal, or any measurement which depends on 
time. But we will also allow functions on lattices such as Z'^. 

In the case = 2, we may be considering the grayscale num- 
bers which represent exposure in a digital camera. In this case, 
the function (grayscale) is defined on a subset of Z^, and is then 
simply a matrix. 

A random walk on Z"* is an assignment of a sequential and 
random motion as a function of time. The randomness presup- 
poses assigned probabilities. But we will use the term "random 
walk" also in connection with random walks on combinatorial 
trees. 

nested refinement multiresolution scales of visual 

subspaces resolutions 

While finite or infinite families of nested subspaces are ubiquitous 
in mathematics, and have been popular in Hilbert space theory 
for generations (at least since the 1930s), this idea was revived in 
a different guise in 1986 by Stephane Mallat, then an engineering 
graduate student. In its adaptation to wavelets, the idea is now 
referred to as the multiresolution method. 

What made the idea especially popular in the wavelet com- 
munity was that it offered a skeleton on which various discrete 
algorithms in applied mathematics could be attached and turned 
into wavelet constructions in harmonic analysis. In fact what we 
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now call multiresolutions have come to signify a crucial link be- 
tween the world of discrete wavelet algorithms, which are popular 
in computational mathematics and in engineering (signal/image 
processing, data mining, etc.) on the one side, and on the other 
side continuous wavelet bases in function spaces, especially in 
L^(R''). Further, the multiresolution idea closely mimics how 
fractals are analyzed with the use of finite function systems. 

But in mathematics, or more precisely in operator theory, 
the underlying idea dates back to work of John von Neumann, 
Norbert Wiener, and Herman Wold, where nested and closed sub- 
spaces in Hilbert space were used extensively in an axiomatic ap- 
proach to stationary processes, especially for time series. Wold 
proved that any (stationary) time series can be decomposed into 
two different parts: The first (deterministic) part can be exactly 
described by a linear combination of its own past, while the sec- 
ond part is the opposite extreme; it is unitary, in the language of 
von Neumann. 

von Neumann's version of the same theorem is a pillar in 
operator theory. It states that every isometry in a Hilbert space 
H is the unique sum of a shift isometry and a unitary operator, 
i.e., the initial Hilbert space Ti. splits canonically as an orthogonal 
sum of two subspaces Hs and Hu in H, one which carries the shift 
operator, and the other W„ the unitary part. The shift isometry 
is defined from a nested scale of closed spaces Vn, such that the 
intersection of these spaces is Hu- Specifically, 

• • • c F_i c Vo c V^i c ^2 c • • • c y„ c y„+i c • • • 
f\Vn=nu, and \JVn=H. 

n n 

However, Stephane Mallat was motivated instead by the no- 
tion of scales of resolutions in the sense of optics. This in turn 
is based on a certain "artificial-intelligence" approach to vision 
and optics, developed earlier by David Marr at MIT, an approach 
which imitates the mechanism of vision in the human eye. 

The connection from these developments in the 1980s back 
to von Neumann is this: Each of the closed subspaces Vn corre- 
sponds to a level of resolution in such a way that a larger subspace 
represents a finer resolution. Resolutions arc relative, not abso- 
lute! In this view, the relative complement of the smaller (or 
coarser) subspace in larger space then represents the visual detail 
which is added in passing from a blurred image to a finer one, i.e., 
to a finer visual resolution. 

This view became an instant hit in the wavelet commu- 
nity, as it offered a repository for the fundamental father and 
the mother functions, also called the scaling function ip, and the 
wavelet function tp. Via a system of translation and scaling op- 
erators, these functions then generate nested subspaces, and we 
recover the scaling identities which initialize the appropriate 
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MATHEMATICS PROBABILITY ENGINEERING PHYSICS 

algorithms. What results is now called the family of pyramid al- 
gorithms in wavelet analysis. The approach itself is called the 
multiresolution approach (MRA) to wavelets. And in the mean- 
time various generalizations (GMRAs) have emerged. 

In all of this, there was a second "accident" at play: As it turned 
out, pyramid algorithms in wavelet analysis now lend themselves 
via multircsolutions, or nested scales of closed subspaces, to an 
analysis based on frequency bands. Here we refer to bands of 
frequencies as they have already been used for a long time in 
signal processing. 

One reason for the success in varied disciplines of the same 
geometric idea is perhaps that it is closely modeled on how we 
historically have represented numbers in the positional number 
system. Analogies to the Euclidean algorithm seem especially 
compelling. 

operator process black box observable 

(if selfadjoint) 

In linear algebra students arc familiar with the distinctions be- 
tween (linear) transformations T (here called "operators") and 
matrices. For a fixed operator T: F — > W, there is a variety of 
matrices, one for each choice of basis in V and in W. In many 
engineering applications, the transformations are not restricted 
to be linear, but instead represent some experiment ( "black box," 
in Norbert Wiener's terminology) , one with an input and an out- 
put, usually functions of time. The input could be an external 
voltage function, the black box an electric circuit, and the output 
the resulting voltage in the circuit. (The output is a solution to 
a differential equation.) 

This context is somewhat different from that of quantum me- 
chanical (QM) operators T: V ^ V where F is a Hilbert space. 
In QM, selfadjoint operators represent observables such as posi- 
tion Q and momentum P, or time and energy. 

Fourier dual generating time/frequency P/Q 

pair function 

The following dual pairs position Q/momentum P, and time/energy 
may be computed with the use of Fourier series or Fourier trans- 
forms; and in this sense they arc examples of Fourier dual pairs. 
If for example time is discrete, then frequency may be represented 
by numbers in the interval [ 0, 27r); or in [ 0, 1) if we enter the num- 
ber 27r into the Fourier exponential. Functions of the frequency 
are then periodic, so the two endpoints are identified. In the case 
of the interval [ 0, 1), on the left is identified with 1 on the right. 
So a low frequency band is an interval centered at 0, while a high 
frequency band is an interval centered at 1/2. 
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MATHEMATICS PROBABILITY ENGINEERING PHYSICS 

Let a function W on [0,1) represent a probability assignment. 
Such functions W are thought of as "filters" in signal process- 
ing. Wc say that W is low-pass if it is 1 at 0, or if it is near 
1 for frequencies near 0. Low-pass filters pass signals with low 
frequencies, and block the others. 

If instead some filter is 1 at 1/2, or takes values near 1 for 
frequencies near 1/2, then we say that W is high-pass; it passes 
signals with high frequency. 



convolution — filter smearing 

Pointwisc multiplication of functions of frequencies corresponds 
in the Fourier dual time-domain to the operation of convolution 
(or of Cauchy product if the time-scale is discrete.) The process 
of modifying a signal with a fixed convolution is called a linear 
filter in signal processing. The corresponding Fourier dual fre- 
quency function is then referred to as "frequency response" or 
the "frequency response function." 

More generally, in the continuous case, since convolution 
tends to improve smoothness of functions, physicists call it "smear- 
ing." 



decomposition — analysis frequency 

(e.g., Fourier components 
coefficients in a 
Fourier expansion) 

Calculating the Fourier coefScients is "analysis," and adding up 
the pure frequencies (i.e., summing the Fourier series) is called 
synthesis. But this view carries over more generally to engineer- 
ing where there are more operations involved on the two sides, 
e.g., breaking up a signal into its frequency bands, transforming 
further, and then adding up the "banded" functions in the end. 
If the signal out is the same as the signal in, we say that the 
analysis/synthesis yields perfect reconstruction. 

integrate reconstruct synthesis superposition 

(e.g., inverse 
Fourier transform) 

Here the terms related to "synthesis" refer to the second half 
of the kind of signal-processing design outlined in the previous 
paragraph. 
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subspace 



resolution 



(signals in a) 
frequency band 

For a space of functions (signals) , the selection of certain frequen- 
cies serves as a way of selecting special signals. When the process 
of scaling is introduced into optics of a digital camera, we note 
that a nested family of subspaces corresponds to a grading of 
visual resolutions. 



Cuntz relations — 



perfect 

reconstruction 
from subbands 



subband 
decomposition 



inner product 



N-l 

SiS* = 1, and S*Sj — 6ijl. 

i=0 

correlation treinsition 
probability 



probability 
of transition 
from one state 
to another 



In many applications, a vector space with inner product captures 
perfectly the geometric and probabilistic features of the situation. 
This can be axiomatized in the language of Hilbert space; and the 
inner product is the most crucial ingredient in the familiar axiom 
system for Hilbert space. 



/out = Tfin — input / output transformation 

of states 

Systems theory language for operators T: F — » W. Then vectors 
in V are input, and in the range of T output. 



fractal — — — 

Intuitively, think of a fractal as reflecting similarity of scales such 
as is seen in fern-like images that look "roughly" the same at 
small and at large scales. Fractals are produced from an infinite 
iteration of a finite set of maps, and this algorithm is perfectly 
suited to the kind of subdivision which is a cornerstone of the 
discrete wavelet algorithm. Self-similarity could refer alternately 
to space, and to time. And further versatility is added, in that 
flexibility is allowed into the definition of "similar." 
— — data mining — 

The problem of how to handle and make use of large volumes of 

data is a corollary of the digital revolution. As a result, the sub- 
ject of data mining itself changes rapidly. Digitized information 
(data) is now easy to capture automatically and to store electron- 
ically. In science, in commerce, and in industry, data represents 
collected observations and information: In business, there is data 
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on markets, competitors, and customers. In manufacturing, there 
is data for optimizing production opportunities, and for improv- 
ing processes. A tremendous potential for data mining exists in 
medicine, genetics, and energy. But raw data is not always di- 
rectly usable, as is evident by inspection. A key to advances is 
our ability to extract information and knowledge from the data 
(hence "data mining"), and to understand the phenomena govern- 
ing data sources. Data mining is now taught in a variety of forms 
in engineering departments, as well as in statistics and computer 
science departments. 

One of the structures often hidden in data sets is some degree 
of scale. The goal is to detect and identify one or more natural 
global and local scales in the data. Once this is done, it is often 
possible to detect associated similarities of scale, much like the fa- 
miliar scale-similarity from multidimensional wavelets, and from 
fractals. Indeed, various adaptations of wavelet-like algorithms 
have been shown to be useful. These algorithms themselves are 
useful in detecting scale-similarities, and are applicable to other 
types of pattern recognition. Hence, in this context, generalized 
multiresolutions offer another tool for discovering structures in 
large data sets, such as those stored in the resources of the Inter- 
net. Because of the sheer volume of data involved, a strictly man- 
ual analysis is out of the question. Instead, sophisticated query 
processors based on statistical and mathematical techniques are 
used in generating insights and extracting conclusions from data 
sets. 



Multiresolutions. Haar's work in 1909 1910 had implicitly the key idea which got 
wavelet mathematics started on a roll 75 years later with Yves Meyer, Ingrid 
Daubechies, Stephane Mallat, and others — namely the idea of a multiresolution. 
In that respect Haar was ahead of his time. See Figures 1 and 2 for details. 

The word "multiresolution" suggests a connection to optics from physics. So 
that should have been a hint to mathematicians to take a closer look at trends 
in signal and image processing! Moreover, even staying within mathematics, it 
turns out that as a general notion this same idea of a "multiresolution" has long 
roots in mathematics, even in such modern and pure areas as operator theory and 
Hilbcrt-spacc geometry. Looking even closer at these interconnections, we can now 
recognize scales of subspaces (so-called multiresolutions) in classical algorithmic 
construction of orthogonal bases in inner-product spaces, now taught in lots of 
mathematics courses under the name of the Gram Schmidt algorithm. Indeed, a 
closer look at good old Gram-Schmidt reveals that it is a matrix algorithm. Hence 
new mathematical tools involving non-commutativity! 

If the signal to be analyzed is an image, then why not select a fixed but suit- 
able resolution (or a subspace of signals corresponding to a selected resolution), 
and then do the computations there? The selection of a fixed "resolution" is dic- 
tated by practical concerns. That idea was key in turning computation of wavelet 
coefficients into iterated matrix algorithms. As the matrix operations get large. 
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the computation is carried out in a variety of paths arising from big matrix prod- 
ucts. The dichotomy, continuous vs. discrete, is quite famihar to engineers. The 
industrial engineers typically work with huge volumes of numbers. 

Numbers! — So why wavelets? Well, what matters to the industrial engineer 
is not really the wavelets, but the fact that special wavelet functions serve as an 
efficient way to encode large data sets — I mean encode for computations. And the 
wavelet algorithms are computational. They work on numbers. Encoding numbers 
into pictures, images, or graphs of functions comes later, perhaps at the very end of 
the computation. But without the graphics, I doubt that we would understand any 
of this half as well as we do now. The same can be said for the many issues that 
relate to the crucial mathematical concept of self-similarity, as we know it from 
fractals, and more generally from recursive algorithms. 

1. Definition 

In this paper we outline several points of view on the interplay between discrete 
and continuous wavelet transforms; stressing both pure and applied aspects of both. 
We outline some new links between the two transform technologies based on the 
theory of representations of generators and relations. By this we mean a finite 
system of generators which are represented by operators in Hilbert space. We 
further outline how these representations yield sub-band filter banks for signal and 
image processing algorithms. 

The word "wavelet transform" (WT) means different things to different people: 
Pure and applied mathematicians typically give different answers the questions 
"What is the WT?" And engineers in turn have their own preferred quite different 
approach to WTs. Still there are two main trends in how WTs are used, the 
continuous WT on one side, and the discrete WT on the other. Here we offer a 
userfriendly outline of both, but with a slant toward geometric methods from the 
theory of operators in Hilbert space. 

Our paper is organized as follows: For the benefit of diverse reader groups, we 
begin with Glossary (section ). This is a substantial part of our account, and it 
reflects the multiplicity of how the subject is used. 

The concept of multiresolutions or multiresolution analysis (MRA) serves as a 
link between the discrete and continuous theory. 

In section |4l we summarize how different mathematicians and scientists have 
contributed to and shaped the subject over the years. 

The next two sections then offer a technical overview of both discrete and the 
continuous WTs. This includes basic tools from Fourier analysis and from operators 
in Hilbert space. In sections[n]and[7]we outline the connections between the separate 
parts of mathematics and their applications to WTs. 

2. Introduction 

While applied problems such as time series, signals and processing of digital 
images come from engineering and from the sciences, they have in the past two 
decades taken a life of their own as an exciting new area of applied mathemat- 
ics. While searches in Google on these keywords typically yield sites numbered in 
the millions, the diversity of applications is wide, and it seems reasonable here to 
narrow our focus to some of the approaches that are both more mathematical and 
more recent. For references, see for example [Tl [6l [23l I31j . In addition, our own 
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interests (e.g., pDl[^[7ni28j ) have colored the presentation below. Each of the two 
areas, the discrete side, and the continuous theory is huge as measured by recent 
journal publications. A leading theme in our article is the independent interest in 
a multitude of interconnections between the discrete algorithm and their uses in 
the more mathematical analysis of function spaces (continuous wavelet transforms) . 
The mathematics involved in the study and the applications of this interaction we 
feel is of benefit to both mathematicians and to engineers. See also .20 . An early 
paper [9] by Daubechies and Lagarias was especially influential in connecting the 
two worlds, discrete and continuous. 



3.1. The Discrete Wavelet Transform. If one stays with function spaces, it is 
then popular to pick the d-dimensional Lebesgue measure on M'*, d = 1,2,, and 
pass to the Hilbert space L^(R'') of all square integrable functions on R'', referring 
to d-dimensional Lebesgue measure. A wavelet basis refers to a family of basis 
functions for L^(IR'^) generated from a finite set of normalized functions V'l , the 
index i chosen from a fixed and finite index set /, and from two operations, one 
called scaling, and the other translation. The scaling is typically specified by a d 
by d matrix over the integers Z such that all the eigenvalues in modulus are bigger 
than one, lie outside the closed unit disk in the complex plane. The d-lattice is 
denoted Z"* , and the translations will be by vectors selected from Z''. We say that 
we have a wavelet basis if the triple indexed family V'ij,fc(a;) := \detA\^/'^'il}{A^ x + k) 
forms an orthonormal basis (ONB) for L^(R'') as i varies in /, j € Z, and k € M''. 
The word "orthonormal" for a family F of vectors in a Hilbert space Ti. refers to 
the norm and the inner product in Ti.: The vectors in an orthonormal family F are 
assumed to have norm one, and to be mutually orthogonal. If the family is also 
total (i.e., the vectors in F span a subspace which is dense in H), we say that F is 
an orthonormal basis (ONB.) 

While there are other popular wavelet bases, for example frame bases, and dual 
bases (see e.g., [21 [T^ and the papers cited there), the ONBs are the most agreeable 
at least from the mathematical point of view. 

That there are bases of this kind is not at all clear, and the subject of wavelets in 
this continuous context has gained much from its connections to the discrete world 
of signal- and image processing. 

Here we shall outline some of these connections with an emphasis on the math- 
ematical context. So we will be stressing the theory of Hilbert space, and bounded 
linear operators acting in Hilbert space Ti., both individual operators, and families 
of operators which form algebras. 

As was noticed recently the operators which specify particular subband algo- 
rithms from the discrete world of signal- processing turn out to satisfy relations 
that were found (or rediscovered independently) in the theory of operator algebras, 
and which go under the name of Cuntz algebras, denoted On if is the number of 
bands. For additional details, see e.g., 

In symbols the C*— algebra has generators (Si)^^^, and the relations are 



3. The discrete vs continuous wavelet Algorithms 



N-l 



(3.1) 
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(where 1 is the identity element in On) and 



N-l 



(3.2) J2 ^^^^ = 1' and S*Sj = (5,^1. 



In a representation on a Hilbert space, say H, the symbols Si turn into bounded 
operators, also denoted Si, and the identity element 1 turns into the identity oper- 
ator / in 7i, i.e., the operator I : h —i- h, for h € Ti.. In operator language, the two 
formulas 13. II and 13.21 state that each Si is an isometry in Ti, and that te respective 
ranges SiTi. are mutually orthogonal, i.e., SiTC _L SjTi. for i ^ j. Introducing the 
projections Pi = SiS* , we get PiPj = 6ijPi, and 



In the engineering literature this takes the form of programming diagrams: 

If the process of Figure 3 is repeated, we arrive at the discrete wavelet transform 
or stated in the form of images (n = 5) 

Selecting a resolution subspace Vq = closure span{ip{- — k)\k G Z}, wc arrive 
at a wavelet subdivision {4'j,k\j > 0,k £ Z}, where 4'j,k{x) — 2^/'^ip{2^x — k), and 
the continuous expansion / = '^jk{'4'j,k\f)4'j.k or the discrete analogue derived 
from the isometrics, i = 1, 2, • • • , TV — 1, S^Si for /c = 0, 1, 2, • • • ; called the discrete 
wavelet transform. 

Notational convention. In algorithms, the letter N is popular, and often used for 
counting more than one thing. 

In the present contest of the Discete Wavelet Algorithm (DWA) or DWT, we 
count two things, "the number of times a picture is decomposed via subdivision". 
We have used n for this. The other related but different number N is the number of 
subbands, = 2 for the dyadic DWT, and A = 4 for the image DWT. The image- 
processing WT in our present context is the tensor product of the 1-D dyadic 
WT, so 2 X 2 = 4. Caution: Not all DWAs arise as tensor products of A^ = 2 
models. The wavelets coming from tensor products are called separable. When a 
particular image-processing scheme is used for generating continuous wavelets it is 
not transparent if we are looking at a separable or inseparable wavelet! 

To clarify the distinction, it is helpful to look at the representations of the Cuntz 
relations by operators in Hilbert space. We are dealing with representations of the 
two distinct algebras ©2, and O4; two frequency subbands vs 4 subbands. Note 
that the Cuntz O2, and ©4 are given axiomatic, or purely symbolically. It is only 
when subband filters are chosen that we get representations. This also means that 
the choice of A^ is made initially; and the same A^ is used in different runs of the 
programs. In contrast, the number of times a picture is decomposed varies from 
one experiment to the next! 

Summary: A^ = 2 for the dyadic DWT: The operators in the representation 
are Sq , Si. One average operator, and one detail operator. The detail operator Si 
"counts" local detail variations. 

Image-processing. Then A^ = 4 is fixed as we run different images in the DWT: 
The operators are now: Sq , Sh, Sy, So- One average operator, and three detail 
operator for local detail variations in the three directions in the plane. 



N-l 
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3.2. The Continuous Wavelet Transform. Consider functions / on the real 
line M. We select the Hilbert space of functions to be L^(M) To start a continuous 
WT, we must select a function -if) e L^(R) and r, s e M. such that the following 
family of functions 

r 

creates an over-complete basis for L^(E). An over-complete family of vectors in 
a Hilbert space is often called a coherent decomposition. This terminology comes 
from quantum optics. What is needed for a continuous WT in the simplest case is 
the following representation valid for all f G (M) : 



where := \ip{Lu)\'^^ and where {ipr,s\f) = J^'>IJr,s{y)fiy)dy. The refinements 
and implications of this are spelled out in tables in section [331 



3.3. Some background on Hilbert space. Wavelet theory is the art of finding a 
special kind of basis in Hilbert space. Let H he a. Hilbert space over C and denote 
the inner product (• | •). For us, it is assumed linear in the second variable. If 
H^L^ (M), then 

If W = £2 (Z), then 



;/|g):= / fix)g{x) 
Jr 



dx. 



Let T = M/27rZ. If 7^ = (t), then 



/(e) 9 (9) dS^ 



inO 



Functions f ^ L'^ (T) have Fourier series: Setting Cn (9) = e 

fin) :=(e„|/) = ^ T ^ (0) df?, 



and 



Il2(t) — / 



Similarly if / G (M), then 

fit) := / e-'-'fix) dx, 

and 

ll/lli^(R) = ^/j/W" 

Let J be an index set. We shall only need to consider the case when J is 
countable. Let {V'alaGj be a family of nonzero vectors in a Hilbert space H. We 
say it is an orthonormal basis (ONB) if 

(3.3) {^pa I "0/3) = Sa,f3 (Kroncckcr delta) 



dt. 
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and if 

(3.4) = 11/11' holds for all / e 

If only (|3.4p is assumed, but not (|3.3p . we say that {ipa\aeJ ^ (normalized) tight 
frame. We say that it is a frame with frame constants 0<yl<_B<ooif 

All/ll' < ^K^a I / <S||/||' holds for aU / e 7i. 

Introducing the rank-one operators Qa '■— \'ipa) (V'al of Dirac's terminology, see [3], 
we see that {^/'aj^gj is an ONB if and only if the Qq's are projections, and 

(3.5) Qa = I (= the identity operator in li). 

aeJ 

It is a (normalized) tight frame if and only if (13. 5p holds but with no further 
restriction on the rank-one operators Qa- It is a frame with frame constants A and 
B if the operator 

satisfies 

AI <S <BI 

in the order of hermitian operators. (We say that operators Hi = H*, i = 1,2, 
satisfy Hi< H2 if if \ Hif ) < ( / | -^2/ ) holds for aU f eH). If ft, are vectors 
in a Hilbert space H, then the operator A = \h) {k\ is defined by the identity 
{u \ Av) ^ {u \ h) {k \ v) foi all u,v e H. 

Wavelets in (M) are generated by simple operations on one or more functions 
ip in (M), the operations come in pairs, say scaling and translation, or phase- 
modulation and translations. If e {2, 3, . . . } we set 

i^j.k (x) := N^l'^il) {N^x - k) for j, k e Z. 

3.3.1. Increasing the dimension. In wavelet theory, there is a tradition for re- 
serving If for the father function and for the mother function. A 1-level wavelet 
transform of an x AI image can be represented as 

, ii I hi' 

(3.6) fi 

where the subimages h^,d^,a^ and each have the dimension of N/2 by Af/2. 



(3.7) 



ai = y„i : ^^{x, y) = ^{x)v{y) = Y.i h,h,^{2x - i)v{2y - j) 

h^=V^® : ^P"{x, y) = ^(xMy) = 9^hM^^ " ^M^V " j) 

^' = W^® : ^P^'ix, y) = <f{x)^{y) = Ej Kg,v{2x - i)^(2y - j) 

d^=W^,® Wl : ^^(x, y) = i,(x)i,{y) = Ej 3.5j¥'(2a; - i)^(2y - j) 

where if is the father function and ip is the mother function in sense of wavelet, 
V space denotes the average space and the W spaces are the difference space from 
multiresolution analysis (MRA) [7]. 

In the formulas, we have the following two indexed number systems a := (hi) 
and d := (gi), a is for averages, and d is for local differences. They are really the 
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input for the DWT. But they also are the key Hnk between the two transforms, the 
discrete and continuous. The hnk is made up of the foUowing scaling identities: 

ip{x) = 2 hi(p{2x — i); 

i'ix) ^ 2^.gj(/7(2x - i); 

and (low-pass normalization) X^iez — ^- The scalars (hi) may be real or complex; 
they may be finite or infinite in number. If there are four of them, it is called the 
"four tap" , etc. The finite case is best for computations since it corresponds to 
compactly supported functions. This means that the two functions ip and ■0 will 
vanish outside some finite interval on a real line. 

The two number systems are further subjected to orthgonality relations, of which 

(3.8) ^ hihi+2k = ^<5o,fc 

is the best known. 

The systems h and g arc both low-pass and high-pass filter coefficients. In 13.71 

denotes the first averaged image, which consists of average intensity values of 
the original image. Note that only ip function, V space and h coefficients are used 
here. Similarly, denotes the first detail image of horizontal components, which 
consists of intensity difference along the vertical axis of the original image. Note 
that if function is used on y and ip function on x, W space for x values and V 
space for y values; and both h and g coefficients are used accordingly. The data 

denotes the first detail image of vertical components, which consists of intensity 
difference along the horizontal axis of the original image. Note that ip function is 
used on x and ip function on y, W space for y values and V space for x values; 
and both h and g coefficients are used accordingly. Finally, denotes the first 
detail image of diagonal components, which consists of intensity difference along 
the diagonal axis of the original image. The original image is reconstructed from 
the decomposed image by taking the sum of the averaged image and the detail 
images and scaling by a scaling factor. It could be noted that only tp function, W 
space and g coefficients are used here. See [33l [28]. 

This decomposition not only limits to one step but it can be done again and 
again on the averaged detail depending on the size of the image. Once it stops at 
certain level, quantization (see [261 132p is done on the image. This quantization 
step may be lossy or lossless. Then the lossless entropy encoding is done on the 
decomposed and quantized image. 

The relevance of the system of identities ([3.8p may be summarized as follows. 
Set 

TOo(z) - hkz'^ for all z G T; 

gk ■■= {-l)''hi-k for aU k <E Z; 

miiz) := ^^gkz''; and 

fcez 

{S,f){z) ^ V2m,iz)f{z^), for j - 0, 1, / G L^{T), z e T. 
Then the following conditions are equivalent: 
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(a) The system of equations (|3.8p is satisfied. 

(b) The operators 5*0 and Si satisfy the Cuntz relations. 

(c) We have perfect reconstruction in the subband system of Figure 3. 

Note that the two operators 5*0 and have equivalent matrix representations. 
Recall that by Parseval's formula we have L^(T) ~ P(Z). So representing So instead 
as an oo X oo matrix acting on column vectors x — {xj)j^z we get 

{Sox)^ = hi^2jXj 
and for the adjoint operator i^o *S'q, we get the matrix representation 

with the overbar signifying complex conjugation. This is computational significance 
to the two matrix representations, both the matrix for ^o, and for Fq := Sq, is 
slanted. However, the slanting of one is the mirror- image of the other, i.e.. 

Significance of slanting. The slanted matrix representations refers to the corre- 
sponding operators in L^. In general operators in Hilbert function spaces have 
many matrix representations, one for each orthonormal basis (ONB), but here we 
are concerned with the ONB consisting of the Fourier frequencies , j d Z. So in 
our matrix representations for the S operators and their adjoints we will be acting 
on column vectors, each infinite column representing a vector in the sequence space 
P. A vector in is said to be of finite size if it has only a finite set of non-zero 
entries. 

It is the matrix Fq that is effective for iterated matrix computation. Reason: 
When a column vector a; of a fixed size, say 2 s is multiplied, or acted on by Fq, 
the result is a vector y of half the size, i.e., of size s. So y ^ Fqx. If we use _Fo Stiid 
Fi together on x, then we get two vectors, each of size s, the other one z = Fix, 
and we can form the combined column vector of y and z] stacking y on top of z. In 
our application, y represents averages, while z represents local differences: Hence 
the wavelet algorithm. 









yo 






yi 




'Fo' 




X-2 
X-l 










xo 

Xl 


Z-l 




X2 


zo 






Zl 







y = Fqx 
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z = F\x 



3.4. Connections to group theory. The first Hne in the two tables below is the 

continuous wavelet transform. It comes from what in physics is called coherent 
vector decompositions. Both transforms applies to vectors in Hilbert space Ti, and 
H may vary from case to case. Common to all transforms is vector input and 
output. If the input agrcH^s with output we say that the combined process yields 
the identity operator image. 1 : H ^ H or written 1-^. So for example if (Si)^^ 
is a finite operator system, and input/output operator example may take the form 

N-l 

SiS* = In- 

i=0 

Summary of and variations on the resolution of the identity operator 1 in or 
in £'^, for ^ and ip where tjjr^s {x) = r~2V (^7^)> 

Cv- = / FT IV' ('^)l^ < oo> 
Jr \^\ 

similarly for and C^^^ = (w) ip (w): 



N = 2 


Ovcrcomplete Basis 


Dual Bases 


continuous 
resolution 








= 1 


= 1 


discrete res- 
olution 


YY\'^j,k) {i^j,k\ = 1, 

j"ez feez 

xpj^k corresponding to 
r = 2-^ s = k2-i 


YY\^j,k) {i>j,k\ =1 

jGZ fcez 


N>2 


Isometrics in £^ 


Dual Operator System in i"^ 


sequence 
spaces 


i=0 

where Sq, ■ ■ ■ , S^-i 
are adjoints to the 
quadrature mirror filter 
operators Fj, i.e., Si = F* 


i=0 

for a dual 
operator system 
Sq, . . . , S'at-i, 
So, . . . , Sn-1 
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Then the assertions in the first table amount to: 





s\f)\' 




) ( 


V'r.s 1 g) 


= ml 


yfeL^ (M) 




(/I.9) 


Vf,geL' (R) 






EE(/i^^-.'^-) 

jez fcGZ 


V^i.fc 1 (, 




= ii/iii. 


VfeL^ (M) 




(/l.9> 




i=Q 




E (^^rci^rd) = 

i=0 







A function ■0 satisfying the resolution identity is called a coherent vector in 
mathematical physics. The representation theory for the (ax + 6)-group, i.e., the 
matrix group G = {(oi) \ ^ b CzM.}, serves as its underpinning. Then the 
tables above illustrate how the {ipj,k} wavelet system arises from a discretization 
of the following unitary representation of G: 

(f,g;,/)W^<.-/(i^) 

acting on (M). This unitary representation also explains 

the discretization step in passing from the first line to the second in the tables 
above. The functions { ipj^k \ j,k £ Z } which make up a wavelet system result from 
the choice of a suitable coherent vector -0 e £^ (M), and then setting 

V'j.fe (x) = (f^(2^^ ''■f')'^') " ^^"^ ^^^^ ~ ■ 

Even though this representation lies at the historical origin of the subject of wave- 
lets, the (ax + 6)-group seems to be now largely forgotten in the next generation 
of the wavelet community. But Chapters 1-3 of [7] still serve as a beautiful pre- 
sentation of this (now much ignored) side of the subject. It also serves as a link to 
mathematical physics and to classical analysis . 



4. List of names and discoveries 
Many of the main discoveries summarized below are now lore. 

1807 Expressing functions as sums of sine and cosine 

Jean Baptiste Joseph waves of frequencies in arithmetic progession (now 

Fourier called Fourier series), 
mathematics, physics 
(heat conduction) 
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1909 

Alfred Haar 
mathematics 



1946 

Denes Gabor 

(Nobel Prize): physics 

(optics, holography) 

1948 

Claude Elwood Shannon 
mathematics, engineering 
(information theory) 

1976 

Claude Garland, Daniel 
Estcban (both) 
signal processing 

1981 

Jean Morlct 
petroleum engineer 



1985 

Yves Meyer 

mathematics, 

applications 

1989 

Albert Cohen 
mathematics (ortho- 
gonality relations), 
numerical analysis 



Discovered, while a student of David Hilbert, an or- 
thonormal basis consisting of step functions, appli- 
cable both to functions on an interval, and functions 
on the whole real line. While it was not realized 
at the time, Haar's construction was a precursor of 
what is now known as the Mallat subdivision, and 
multiresolution method, as well as the subdivision 
wavelet algorithms. 

Discovered basis expansions for what might now be 
called time-frequency wavelets, as opposed to time- 
scale wavelets. 

A rigorous formula used by the phone company for 
sampling speech signals. Quantizing information, 
entropy, founder of what is now called the mathe- 
matical theory of communication. 

Discovered subband coding of digital transmission 
of speech signals over the telephone. 



Suggested the term "ondelettes." J.M. decomposed 
reflected seismic signals into sums of "wavelets (Fr.: 
ondelettes) of constant shape," i.e., a decomposi- 
tion of signals into wavelet shapes, selected from a 
library of such shapes (now called wavelet series). 
Received somewhat late recognition for his work. 
Due to contributions by A. Grossman and Y. Meyer, 
Morlet's discoveries have now come to play a central 
role in the theory. 

Mentor for A. Cohen, S. Mallat, and other of the 
wavelet pioneers, Y.M. discovered infinitely often 
differentiable wavelets. 



Discovered the use of wavelet filters in the analy- 
sis of wavelets — the so-called Cohen condition for 
orthogonality. 
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1986 

Stephane Mallat 
mathematics, signal 
and image processing 



Discovered what is now known as the subdivision, 
and muhiresolution method, as weU as the subdivi- 
sion wavelet algorithms. This allowed the effective 
use of operators in the Hilbert space and of 

the parallel computational use of recursive matrix 
algorithms. 



1987 

Ingrid Daubechies 
mathematics, physics, 
and communications 



Discovered differentiable wavelets, with the number 
of derivatives roughly half the length of the sup- 
port interval. Further found polynomial algorithmic 
for their construction (with coauthor Jeff Lagarias; 
joint spectral radius formulas). 



1991 

Wayne Lawton 
mathematics 
(the wavelet 
transfer operator) 



Discovered the use of a transfer operator in the anal- 
ysis of wavelets: orthogonality and smoothness. 



1992 

The FBI 

using wavelet algo- 
rithms in digitizing and 
compressing 
fingerprints 



C. Brislawn and his group at Los Alamos created 
the theory and the codes which allowed the com- 
pression of the enormous FBI fingerprint file, creat- 
ing A/D, a new database of fingerprints. 



2000 

The International 

Standards 

Organization 



A wavelet-based picture compression standard, 
called JPEG 2000, for digital encoding of images. 



1994 

David Donoho 

statistics, 

mathematics 



Pioneered the use of wavelet bases and tools from 
statistics to "denoise" images and signals. 



5. History 

While wavelets as they have appeared in the mathematics literature (e.g., [7]) 
for a long time, starting with Haar in 1909, involve function spaces, the connec- 
tions to a host of discrete problems from engineering is more subtle. Moreover 
the deeper connections between the discrete algorithms and the function spaces of 
mathematical analysis are of a more recent vintage, see e.g., [31] and [2T] . 

Here we begin with the function spaces. This part of wavelet theory refers to 
continous wavelet transforms (details below). It dominated the wavelet literature 
in the 1980s, and is beautifully treated in the first four chapters in and in [8]. 
The word "continuous" refers to the continuum of the real line M. Here we consider 
spaces of functions in one or more real dimensions, i.e., functions on the line M 
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(signals), the plane (images), or in higher dimensions M'', functions of d real 
variables. 

6. Tools from Mathematics 

In our presentation, we will rely on tools from at least three separate areas of 
mathematics, and we will outline how they interact to form a coherent theory, and 
how they come together to form a link between what is now called the discrete 
and the continuous wavelet transform. It is the discrete case that is popular with 
engineers ( [H EH EH ED] ) , while the continuous case has come to play a central role 
in the part of mathematics referred to as harmonic analysis, [8J. The three areas 
are, operator algebras, dynamical systems, and basis constructions: 

a. Operator algebras. The theory of operator algebras in turn breaks up in 
two parts: One the study of "the algebras themselves" as they emerge 
from the axioms of von Neumann (von Neumann algebras) , and Gelfand, 
Kadison and Segal (C*-algebras.) The other has a more applied slant: 
It involves "the representations" of the algebras. By this we refer to the 
following: The algebras will typically be specified by generators and by 
relations, and by a certain norm-completion, in any case by a system of 
axioms. This holds both for the norm-closed algebras, the so called C*- 
algebras, and for the weakly closed algebras, the von Neumann algebras. 
In fact there is a close connection between the two parts of the theory: For 
example, representations of C*-algebras generate von Neumann algebras. 

To talk about representations of a fixed algebra say A we must specify a 
Hilbert space, and a homomorphism p from A into the algebra B{H) of all 
bounded operators on Ti.. We require that p sends the identity element in 
A into the identity operator acting on Ti, and that p{a*) = {p{a))* where 
the last star now refers to the adjoint operator. 

It was realized in the last ten years (see for example O EU E2] that 
a family of representations that wavelets which are basis constructions in 
harmonic analysis, in signal/image analysis, and in computational math- 
ematics may be built up from representations of an especially important 
family of simple C*-algebras, the Cuntz algebras. The Cuntz algebras are 
denoted 02,03,..., including Ooo- 

b. Dynamical systems. The connection between the Cuntz algebras On for 
iV = 2, 3, are relevant to the kind of dynamical systems which are built 
on branching-laws, the case of O^v representing iV-fold branching. The 
reason for this is that if N is fixed. On includes in its definition an iterated 
subdivision, but within the context of Hilbert space. For more details, see 

e.g., [miiaiisiiiiiiiHiiniES]. 

c. Analysis of bases in function spaces. The connection to basis constructions 
using wavelets is this: The context for wavelets is a Hilbert space H, where 
H may be L^(R'') where d is a dimension, d = 1 for the line (signals), d = 2 
for the plane (images), etc. The more successful bases in Hilbert space are 
the orthonormal bases ONBs, but until the mid 1980s, there were no ONBs 
in L^(R'') which were entirely algorithmic and effective for computations. 
One reason for this is that the tools that had been used for 200 years 
since Fourier involved basis functions (Fourier wave functions) which ere 
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not localized. Moreover these existing Fourier tools were not friendly to 
algorithmic computations. 

7. A Transfer Operator 

A popular tool for deciding if a candidate for a wavelet basis is in fact an ONB 
uses a certain transfer operator. Variants of this operator is used in diverse areas 
of applied mathematics. It is an operator which involves a weighted average over 
a finite set of possibilities. Hence it is natural for understanding random walk 
algorithms. As remarked in for example [20l [211 IHl US], it was also studied in 
physics, for example by David Ruelle who used to prove results on phase transition 
for infinite spin systems in quantum statistical mechanics. In fact the transfer 
operator has many incarnations (many of them known as Ruelle operators) , and all 
of them based on A^-fold branching laws. 

In our wavelet application, the Ruelle operator weights in input over the N 
branch possibilities, and the weighting is assigned by a chosen scalar function W. 
the and the H^- Ruelle operator is denoted Rw- In the wavelet setting there is in 
addition a low-pass filter function tuq which in its frequency response formulation 
is a function on the d-torus T'' = R'^/Z'^. 

Since the scaling matrix A has integer entries A passes to the quotient W^/Z"^, 
and the induced transformation va - T^ T'^ is an A''-fold cover, where N — \detA\, 
i.e., for every x in T'' there are N distinct points y in T'' solving rA{y) = x. 

In the wavelet case, the weight function W is W = |mop. Then with this 
choice of W, the ONB problem for a candidate for a wavelet basis in the Hilbert 
space L^(M'^) as it turns out may be decided by the dimension of a distinguished 
eigenspace for Rw, by the so called Perron- Frobenius problem. 

This has worked well for years for the wavelets which have an especially simple 
algorithm, the wavelets that are initialized by a single function, called the scaling 
function. These are called the multiresolution analysis (MRA) wavelets, or for 
short the MRA-wavelets. But there are instances, for example if a problem must be 
localized in frequency domain, when the MRA-wavelets do not suffice, where it will 
by necessity include more than one scaling function. And we are then back to trying 
to decide if the output from the discrete algorithm, and the On representation is 
an ONB, or if it has some stability property which will serve the same purpose, in 
case where asking for an ONB is not feasible. 

8. Future Directions 

The idea of a scientific analysis by subdividing a fixed picture or object into its 
finer parts is not unique to wavelets. It works best for structures with an inherent 
self-similarity; this self-similarity can arise from numerical scaling of distances. But 
there are more subtle non-linear self-similarities. The Julia sets in the complex plane 
are a case in point [H [51 [TUl HH [Ml [IS] . The simplest Julia set come from a one 
parameter family of quadratic polynomials ifidz) = z'^ + c, where z is a complex 
variable and where c is a fixed parameter. The corresponding Julia sets Jc have 
a surprisingly rich structure. A simple way to understand them is the following: 
Consider the two brances of the inverse P± ^ z i-^ ±-\/z — c. Then Jc is the unique 
minimal non-empty compact subset of C, which is invariant under {/3±}. (There 
are alternative ways of presenting Jc but this one fits our purpose. The Julia set 
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J of a holomorphic function, in this case z i-^ + c, informally consists of those 
points whose long-time behavior under repeated iteration , or rather iteration of 
substitutions, can change drastically under arbitrarily small perturbations.) Here 
"long-time" refers to largen n, where (p'-"-^^^{z) — ip{ip'-"-^ (z)) , n — 0,1,..., and 
(^(0)(z) = z. 

It would be interesting to adapt and modify the Haar wavelet, and the other 
wavelet algorithms to the Julia sets. The two papers [ini [TT] initiate such a devel- 
opment. 

9. Literature 

As evidenced by a simple Google check, the mathematical wavelet literature 
is gigantic in size, and the manifold applications spread over a vast number of 
engineering journals. While we cannot do justice to this volumest literature, we 
instead offer a collection of the classics [TH] edited recently by C. Heil et.al. 

Acknowledgements. We thank Professors Dorin Dutkay and Judy Packer for helpful 
discussions. 

10. Bibliography 



[1] Gilles Aubert and Pierre Kornprobst. Mathematical problems in image processing. 2006. 
[2] Lawrence Baggett, Palle Jorgensen, Kathy Merrill, and Judith Packer. A non-MRA C frame 

wavelet with rapid decay. Acta Appl. Math., 2005. 
[3] 0. Bratelli and P. Jorgensen. Wavelets Through a Looking Glass: The World of the Spectrum. 

Birkhiiuser, 2002. 

[4] M. Braverman and M. Yampolsky. Non-computable Julia sets. J. Amer. Math. Soc, 

19(3):551-578 (electronic), 2006. 
[5] Mark Braverman. Parabolic Julia sets are polynomial time computable. Nonlinearity, 

19(6):1383-1401, 2006. 

[6] K. Bredies, D. A. Lorenz, and P. Maass. An optimal control problem in medical image pro- 
cessing. 2006. 

[7] Ingrid Daubechies. Ten lectures on wavelets, volume 61 of CBMS-NSF Regional Conference 

Series in Applied Mathematics. 1992. 
[8] Ingrid Daubechies. Wavelet transforms and orthonormal wavelet bases. Proc. Sympos. Appl. 

Math. 1993. 

[9] Ingrid Daubechies and Jeffrey C. Lagarias. Two-scale difference equations. II. Local regularity, 
infinite products of matrices and fractals. SIAM J. Math. Anal., 1992. 

[10] Robert L. Devaney and Daniel M. Look. A criterion for Sierpinski curve Julia sets. Topology 
Proc, 30(1):163-179, 2006. Spring Topology and Dynamical Systems Conference. 

[11] Robert L. Devaney, Monica Moreno Rocha, and Stefan Siegmund. Rational maps with gen- 
eralized Sierpinski gasket Julia sets. Topology Appl., 154(1):11— 27, 2007. 

[12] Dorin E. Dutkay and Palle E. T. Jorgensen. Wavelets on fractals. Rev. Mat. Iberoamericana, 
22, 2006. 

[13] Dorin E. Dutkay and Kjetil Roysland. The algebra of harmonic functions for a matrix-valued 

transfer operator. arXiv:math/0611539, 2007. 
[14] Dorin E. Dutkay and Kjetil Roysland. Covariant representations for matrix-valued transfer 

operators. arXiv:math/070145.3, 2007. 
[15] Dorin Ervin Dutkay. The spectrum of the wavelet Galerkin operator. Integral Eguations 

Operator Theory, 2004. 

[16] Dorin Ervin Dutkay and Palle E. T. Jorgensen. Wavelet constructions in non-linear dynamics. 

Electron. Res. Announc. Amer. Math. Soc, 2005. 
[17] Dorin Ervin Dutkay and Palle E. T. Jorgensen. Hilbert spaces built on a similarity and on 

dynamical rcnormalization. J. Math. Phys., 2006. 



22 



PALLE E.T. JORGENSEN AND MYUNG-SIN SONG 



[18] Dorin Ervin Dutkay and Palle E. T. Jorgensen. Iterated function systems, Ruelle operators, 
and invariant projective measures. Math. Comp., 2006. 

[19] Christopher Heil and David F. Walnut, editors. Fundamental papers in wavelet theory. Prince- 
ton University Press, Princeton, NJ, 2006. 

[20] Pallc E. T. Jorgensen. Matrix factorizations, algorithms, wavelets. Notices Amer. Math. Soc., 
50, 2003. 

[21] Palle E. T. Jorgensen. Analysis and probability: wavelets, signals, fractals, volume 234 of 

Graduate Texts in Mathematics. Springer, New York, 2006. 
[22] Palle E. T. Jorgensen. Certain representations of the Cuntz relations, and a question on 

wavelets decompositions. 414:165-188, 2006. 
[23] Feng Liu. Diffusion filtering in image processing based on wavelet transform. Sci. China Ser. 

F, 49, 2006. 

[24] John Milnor. Pasting together Julia sets: a worked out example of mating. Experiment. 

Math., 13{l):55-92, 2004. 
[25] C. L. Petersen and S. Zakcri. On the Julia set of a typical quadratic polynomial with a Sicgcl 

disk. Ann. of Math. (2), 159(l):l-52, 2004. 
[26] A. Skodras, C. Christopoulos, and T. Ebrahimi. Jpeg 2000 still image compression standard" 

ieee signal processing magazine. IEEE Signal processing Magazine, 18:36—58, Sept. 2001. 
[27] Myung-Sin Song. Wavelet image compression. Ph.D. thesis. The University of Iowa, 2006. 
[28] Myung-Sin Song. Wavelet image compression. In Operator theory, operator algebras, and 

applications, volume 414 of Contemp. Math., pages 41-73. Amer. Math. Soc, Providence, 

RI, 2006. 

[29] Gilbert Strang. Wavelets from filter banks. Springer, 1997. 

[30] Gilbert Strang. Signal processing for everyone, volume 1739 of Lecture Notes in Math. 2000. 
[31] Gilbert Strang and Truong Nguyen. Wavelets and filter banks. Wellesley-Cambridge Press, 

WcUcslcy, MA, 1996. 

[32] B. E. Uscvitch. A tutorial on modern lossy wavelet image compression: Foundations of jpeg 

2000. IEEE Signal processing Magazine, 18:22-35, Sept. 2001. 
[33] J. S Walker. A Primer on Wavelets and their Scientific Applications. Chapman & Hall, CRC, 

1999. 

(Palle E.T. Jorgensen) Department of Mathematics, The University of Iowa, 14 MacLean 
Hall, Iowa City, IA 52242 

E-mail address: jorgenamath.uiowa.edu 

(Myung-Sin Song) Department of Mathematics and Statistics, Southern Illinois Uni- 
versity Edwardsville, Box 1653, Science Building, Edwardsville, IL 62026 
E-mail address: msongSs iue.edu 



